参考: http://quant-econ.net/py/optgrowth.html
http://www.rieb.kobe-u.ac.jp/academic/ra/dp/English/dp180.pdf
http://www.econ.hokudai.ac.jp/~kudoh/Bellman.pdf
http://www.hara.kier.kyoto-u.ac.jp/NoteByInami3.pdf
Stokey Lucas with Prescott Chapter 2, 3, 4, (9)
In [34]:
%matplotlib inline
import math
import numpy as np
from scipy.optimize import fminbound
from scipy import interp
import matplotlib as mpl
import matplotlib.pyplot as plt
import quantecon as qe
np.set_printoptions(precision=5)
plt.rcParams['font.size'] = 13
# 日本語対応
mpl.rcParams['font.family'] = 'Osaka'
次のように, 今日(t=0)と明日(t=1)に分けてケーキを食べる問題を考える.
$$ \begin{align} \max_{c_0, c_1} \hspace{5pt} & u(c_0) + \beta u(c_1) \\ s.t. \ & c_0 + c_1 \leq \bar{x}\\ & c_0 , c_1 \geq 0 \end{align} $$例えば今, $u(c) = \log(c), \ \beta = 0.5, \ \bar{x} = 1$ と置くと, 問題は
$$ \begin{align} \max_{c_0} \hspace{5pt} \log(c_0) + 0.5 \log(1 - c_0) \hspace{5pt} s.t. \hspace{5pt} c_0 \geq 0 \end{align} $$と書き直せるので(ケーキを残すよりは全て消費した方が効用が上がるので, $c_0 + c_1 = 1$が望ましい. 厳密にはKKT条件を使う), これを解くと,
$$ \begin{align} \it{L} &= \log(c_0) + 0.5 \log(1 - c_0) \\ \\ \text{F.O.C.は, } &\\ \\ \frac{d\it{L}}{dc_0} & = \frac{1}{c_0} - \frac{0.5}{1-c_0} = 0 \\ \therefore c_0 & = \frac{2}{3} \end{align} $$となって, 今日2/3, 明日1/3消費するのが最適になる.
同じ問題をn期間で考える.
$$ \begin{align} \max_{c_t, \ t \in \{0, 1, \ldots, n-1 \}} \hspace{5pt} & \sum_{t=0}^{n-1} \beta^{t} u(c_t) \\ s.t. \ & c_0 + \ldots + c_{n-1} \leq \bar{x}\\ & c_0 , \ldots, c_{n-1} \geq 0 \end{align} $$ラグランジュの未定乗数法で解く.
$$ \begin{align} \it{L} &= \sum_{t=0}^{n-1} \beta^{t} u(c_t) + \lambda (c_0 + \ldots + c_{n-1} - \bar{x}) \\ \\ \text{F.O.C.は, } &\\ \\ \frac{\partial \it{L}}{\partial c_t} & = \beta^t u'(c_t) + \lambda = 0 \hspace{15pt} \forall t \in \{0, 1, \ldots, n-1 \} \\ \frac{\partial \it{L}}{\partial \lambda} & = \sum_{t=0}^{n-1} c_t - \bar{x} = 0 \end{align} $$上と同様に, $u(c) = \log(c), \ \beta = 0.5, \ \bar{x} = 1$ と置くと, 第1式から,
$$ \begin{align} \frac{1}{c_0} &= \frac{\beta}{c_1} = \ldots = \frac{\beta^{n-1}}{c_{n-1}} \\ \therefore c_1 &= \beta c_0, c_2 = \beta^2 c_0, \ldots, c_{n-1} = \beta^{n-1} c_0 \end{align} $$となり,
$$ \begin{align} \sum_{t=0}^{n-1} c_t &= \sum_{t=0}^{n-1} \beta^t c_0 \\ &= \frac{1 - \beta^n}{1 - \beta} c_0 \\ &= \frac{1 - 0.5^n}{0.5} c_0 = 1 \\ \therefore c_0 & = \frac{2}{1 - 0.5^n} \end{align} $$となり, 他の$c_1, \ldots, c_{n-1}$ も求まる.
2の問題を, t期開始時点でのケーキの残量を表す変数 $x_0, x_1, \ldots, x_n $ を導入して書き換えることを考えると,
$$ \begin{align} \max_{c_t, x_{t+1}} \hspace{5pt} & \sum_{t=0}^{n-1} \beta^{t} u(c_t) \\ s.t. \ & c_t + x_{t+1} = x_t \hspace{15pt} \forall t \in \{0, 1, \ldots, n-1 \} \\ & x_0 = \bar{x} \\ & c_0 , c_1, \ldots, c_{n-1}, x_0, x_1, \ldots x_n \geq 0 \end{align} $$となる. ただし, 2で与えていた条件 $\sum_{t=0}^{n-1} c_t = \bar{x}$ は除いている. ここで, $c_t$を制御変数(Control Variable), $x_t$を状態変数(State Variable)という.
$c_t$ を消去すると,
$$ \begin{align} \max_{x_0, \ldots, x_n} \hspace{5pt} & \sum_{t=0}^{n-1} \beta^{t} u(x_t - x_{t+1}) \\ s.t. \ & x_0 = \bar{x} \\ & x_0, x_1, \ldots, x_n \geq 0 \end{align} $$と, 制御変数だけで書き直せる. この問題のラグランジュ関数は,
$$ \begin{align} {\it L} = \sum_{t=0}^{n-1} \beta^{t} u(x_t - x_{t+1}) \end{align} $$となり, F.O.Cは,
$$ \begin{align} \frac{\partial \it{L}}{\partial x_{t+1}} &= - \beta^t u'(x_t - x_{t+1}) + \beta^{t+1} u'(x_{t+1} - x_{t+2}) = 0 \hspace{15pt} \forall t \in \{0, 1, \ldots, n-2 \} \\ \frac{\partial \it{L}}{\partial x_n} &= - \beta^{n-1} u'(x_{n-1} - x_{n}) = 0 \end{align} $$ここで, 2番目の条件は
$$ \begin{align} \beta^{n-1} u'(x_{n-1} - x_{n}) = 0 \end{align} $$または
$$ \begin{align} \beta^{n-1} u'(x_{n-1} - x_{n}) \leq 0, x_n = 0 \end{align} $$の時に成り立つ(KKT conditionによる).
2式の意味は「最後の期には資産(ケーキ)を残さない」 か 「残っていたとしても, 無価値」 である状態が最適ということ(今の設定の場合, 残ったケーキが無価値ということはあり得ないが, 無限期間に拡張した場合にはあり得るのでこうしておく).
2つをまとめて
$$ \begin{align} \beta^{n-1} \frac{\partial u(x_{n-1} - x_{n})}{\partial x_n} x_n = 0 \end{align} $$を横断性条件(Transversality Condition)という.
同じ問題を無限期間で考える.
$$ \begin{align} \max_{c_t, \ t \in \{0, 1, \ldots \}} \hspace{5pt} & \sum_{t=0}^{\infty} \beta^{t} u(c_t) \\ s.t. \ & \sum_{t=0}^\infty c_t \leq \bar{x}\\ & c_0 , c_1, \ldots \geq 0 \end{align} $$3と同様に, t期開始時点でのケーキの残量を表す変数 $x_0, x_1, \ldots$ を導入すると, 問題は
$$ \begin{align} \max_{c_t, x_{t+1}} \hspace{5pt} & \sum_{t=0}^{\infty} \beta^{t} u(c_t) \\ s.t. \ & c_t + x_{t+1} = x_t \hspace{15pt} \forall t \in \{0, 1, \ldots \} \\ & x_0 = \bar{x} \\ & c_0 , c_1, \ldots, x_0, x_1, \ldots \geq 0 \end{align} $$と書き直せる.
この問題のラグランジュ関数は,
$$ \begin{align} \it{L} = \sum_{t=0}^{\infty} \beta^{t} \left( u(c_t) + \lambda_t (c_t + x_{t+1} - x_t) \right) \end{align} $$となり, F.O.Cは,
$$ \begin{align} \frac{\partial \it{L}}{\partial c_t} &= \beta^t \left( u'(c_t) + \lambda_t \right) = 0 \\ \frac{\partial \it{L}}{\partial x_{t+1}} &= \beta^t \lambda_t - \beta^{t+1} \lambda_{t+1} = 0 \\ \frac{\partial \it{L}}{\partial \lambda_t} &= c_t + x_{t+1} - x_t = 0 \\ \end{align} $$となる. 変形して,
$$ \begin{align} u'(c_t) &= - \lambda_t \\ \lambda_t &= \beta \lambda_{t+1} \\ \end{align} $$1式を2式に代入して,
$$ \begin{align} & u'(c_t) = \beta u'(c_{t+1}) \\ \text{または, } \hspace{10pt} &\\ & u'(x_{t} - x_{t+1}) = \beta u'(x_{t+1} - x_{t+2}) \end{align} $$を得る. これがオイラー方程式.
3と同様に, オイラー方程式のみからは最適経路は求まらず, 次の横断性条件(Transversality Condition)
$$ \begin{align} \lim_{t \to \infty} \beta^t \frac{\partial u(x_t, x_{t+1})}{\partial x_{t+1}} x_{t+1} = 0\\ \end{align} $$が必要になる. これはn期間の横断性条件の極限をとったものになっている.
オイラー方程式と横断性条件が与えられたとして, どのようにして具体的な解を求めるか?
そもそも, 無限期間のSequential Problemの解が解析的に書ける(例えば$c_t = log(x_t) + 3$ のように )ことはめったにない.
しかし, ごく限られた問題については, 簡単な関数形で統一的に解が書けることが知られている. このような場合には「Guess and Verify」:
を行うことができる.
同じ問題を考える.
$$ \begin{align} \max_{c_t, x_{t+1}} \hspace{5pt} & \sum_{t=0}^{\infty} \beta^{t} u(c_t) \\ s.t. \ & c_t + x_{t+1} = x_t \hspace{15pt} \forall t \in \{0, 1, \ldots \} \\ & x_0 = \bar{x} \\ & c_0 , c_1, \ldots, x_0, x_1, \ldots \geq 0 \end{align} $$はSequential Problemと呼ばれ, $x_0, x_1, \ldots$ の差分方程式になっている.
一方でこの問題は, 価値関数(Value Function)を用いて, 次のように書き換えることができる.
$$ \begin{align} v(x_t) &= \max_{c_t} \hspace{3pt} u(x_t - x_{t+1}) + \beta v(x_{t+1}) \\ s.t. \ & x_0 = \bar{x}\\ & c_0 , c_1, \ldots, x_0, x_1, \ldots \geq 0 \end{align} $$これは関数方程式(Functional Equation)で, $v(x_t)$ は t期時点から見た, t+1期以降の効用の和の最大値を表している.
ここで, 政策関数(Policy Function)
$$ \begin{align} c_t = \sigma_t(x_t) \end{align} $$を考える. これは, t期のState Variable(例えば, t期開始時のケーキの残量)が与えられたもとで, Control Variable(例えば, t期にケーキをどのくらい消費するか)を決定する関数になっている.
一定の仮定のもとで, 政策関数はtに依存せず
$$ \begin{align} c_t = \sigma(x_t) \end{align} $$と一意に表現できる(SLP Chap4).
このことから, 関数方程式の解は
$$ \begin{align} v(x) &= \max_{c_t} \hspace{3pt} u(c_t) + \beta v(x) \end{align} $$の解として表現できる.
この関数方程式と, 横断性条件
$$ \begin{align} \lim_{t \to \infty} \beta^t v(x_t) = 0\\ \end{align} $$を満たす解v(x)は, 元の関数方程式を満たす(SLP Theorem4.5).
一定の仮定のもとで, Sequential Problemの解はFunctional Equationを満たす(SLP Theorem4.4).
一方で, Functional Equationと横断性条件を満たすfeasible planは, Sequential Problemの解を与える.
では, 何故SP以外にFEが必要なのか?
SPからオイラー方程式を導出するためには, 目的関数や制約条件が微分可能でなければならない. また, 最適性を保証するためにはuやfが凹関数であることも必要.
このような条件が満たされない時には, オイラー方程式と横断性条件による解法は使えない. 一方で, FEは各期の最適化の構造が1期間の差分方程式で書けさえすればよい. その点で, FEはより柔軟に問題を表現できる.
ベルマン方程式を解く具体的な方法としては, 大きく分けて
がある.
Guess and Verifyはオイラー方程式の場合とほぼ同じ. Value FunctionまたはPolicy Functionの具体的な関数形を推測し, 1. 関数方程式を満たすか(つまり, 関数方程式を満たすように推測した関数形の係数を決定できるか) 2. 横断性条件を満たすか をチェックすれば良い.
ただし, これが使える問題は非常に限られている(多くの問題ではOptimal Value FunctionやPolicy Functionを簡単な関数形で表現することはできないので). 一般には, 次のIterationによる方法を用いる.
ここでは, Value Function Iteration という方法を説明する. やり方はいたって簡単.
更新式は, Bellman Operator $T$ を用いて, 次のようにも書くこともある.
$$ \begin{align} v_{t+1} &= T v_t \\ \end{align} $$6で扱う最適成長モデルの場合, $u, f$がcontinuousかつ$u$がboundedなどの条件を満たせば, 少なくとも1つの最適解に収束することが示されている(SLP chap4).
次のような単純な最適成長モデルを考える.
$$ \begin{align} \max_{c_t, k_{t+1}} \hspace{5pt} & \sum_{t=0}^{\infty} \beta^{t} u(c_t) \\ s.t. \ & c_t + k_{t+1} = f(k_t) \hspace{15pt} \forall t \in \{0, 1, \ldots \} \\ & k_0 = \bar{k} \\ & c_0 , c_1, \ldots, k_0, k_1, \ldots \geq 0\\ \end{align} $$このSequential ProblemをFunctional Equationの形で書き直すと,
$$ \begin{align} v(k_t) = & \max_{k_{t+1}} \hspace{3pt} u(f(k_t) - k_{t+1}) + \beta v(k_{t+1}) \\ s.t. \ & k_0 = \bar{k}\\ & c_0 , c_1, \ldots, k_0, k_1, \ldots \geq 0\\ \end{align} $$となる. これを解いて, Optimal Policy = 「各期開始時に資本 $k_t$ が与えられたもとで, その資本を$c_t$: 消費 と $x_{t+1}$: 生産 にどう割り振ればいいかを決める政策」を求めることが目標となる.
今, $u(c_t) = log c_t$, $f(k_t) = k_t^\alpha$ と特定化する. この時, Optimal Value Functionは,
$$ v^*(k_t) = A + B \log k_t, \hspace{15pt} A, B: \text{const} $$になることが知られている. 実際, これをFEに代入すると,
$$ v(k_t) = \max_{k_{t+1}} \hspace{3pt} \log(k_t^\alpha - k_{t+1}) + \beta ( A + B \log k_{t+1} ) \\ $$となり, $k_{t+1}$ で微分した条件を用いて係数比較を行うと,
$$ \begin{align} A &= \frac{\log(1-\alpha \beta)}{1 - \beta} + \frac{\log(\alpha \beta) \alpha \beta }{(1 - \alpha \beta)(1 - \beta)}\\ B &= \frac{\alpha}{1 - \alpha \beta} \end{align} $$となる. これが横断性条件を満たすことも示せる(が, 少し面倒).
In [41]:
# Primitives and grid
alpha = 0.65
beta = 0.95
grid_min = 1e-6
grid_max = 5
grid_size = 150
grid = np.linspace(grid_min, grid_max, grid_size)
# Exact solution
ab = alpha * beta
c1 = (np.log(1 - ab) + np.log(ab) * ab / (1 - ab)) / (1 - beta)
c2 = alpha / (1 - ab)
def v_star(k):
return c1 + c2 * np.log(k)
def bellman_operator(w):
# Apply linear interpolation to w
Aw = lambda x: interp(x, grid, w)
# set Tw[i] equal to max_c { log(c) + beta w(f(k_i) - c)}
Tw = np.empty(grid_size)
for i, k in enumerate(grid):
objective = lambda c: - np.log(c) - beta * Aw(k**alpha - c)
c_star = fminbound(objective, grid_min, k**alpha)
Tw[i] = - objective(c_star)
return Tw
if __name__ == '__main__':
w = 5 * np.log(grid) - 25 # An initial condition -- fairly arbitrary
n = 75
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_ylim(-40, -20)
ax.set_xlim(np.min(grid), np.max(grid))
ax.plot(grid, w, color=plt.cm.jet(0), linewidth=2, label="initial Value Function")
for i in range(n-1):
w = bellman_operator(w)
ax.plot(grid, w, color=plt.cm.jet(i / n), linewidth=1, alpha=0.5)
w = bellman_operator(w)
ax.plot(grid, w, color="red", linewidth=2, label="VF after {} repeats".format(n))
ax.plot(grid, v_star(grid), 'k-', linewidth=2, label="optimal Vaule Function")
ax.legend(loc='upper left')
plt.show()
In [38]:
def compute_greedy(w):
Aw = lambda x: interp(x, grid, w)
sigma = np.empty(grid_size)
for i, k in enumerate(grid):
objective = lambda c: - np.log(c) - beta * Aw(k**alpha - c)
sigma[i] = fminbound(objective, 1e-6, k**alpha)
return sigma
In [40]:
true_sigma = (1 - alpha * beta) * grid**alpha
fig, ax = plt.subplots(3, 1, figsize=(8, 10))
for i, n in enumerate((2, 4, 6)):
ax[i].set_ylim(0, 1)
ax[i].set_xlim(0, 2)
ax[i].set_yticks((0, 1))
ax[i].set_xticks((0, 2))
w = 5 * np.log(grid) - 25 # Initial condition
v_star = qe.compute_fixed_point(bellman_operator, w,
max_iter=n,
verbose=0)
sigma = compute_greedy(v_star)
ax[i].plot(grid, sigma, 'b-', lw=2, alpha=0.8, label='approximate optimal policy')
ax[i].plot(grid, true_sigma, 'k-', lw=2, alpha=0.8, label='true optimal policy')
ax[i].legend(loc='upper left')
ax[i].set_title('{} value function iterations'.format(n))
In [ ]: